<?php

declare(strict_types=1);

namespace DataCube\DataCubeAggregation\Functions\Statistic;

class StatisticalGraphics
{
    protected $values    = array();
    protected $dataset   = array();
    protected $precision = 3;

    public function setPrecision($n)
    {
        $this->precision = (int) $n;
    }

    public function getPrecision($n)
    {
        return $this->precision;
    }

    public function boxplot(array $x)
    {
        sort($x);
        $staticFunc = (new StatisticalFunc());
        $count  = count($x);
        $median = $staticFunc->median($x);
        $sd     = $staticFunc->sd($x);
        $q1     = $x[round($count / 4)];
        $q3     = $x[round($count * 3 / 4)];
        $min    = $x[0];
        $max    = $x[$count - 1];

        // Outliars ]median + 2sd, median - 2sd[
        $outliers = array();

        if ($min < $median - 2 * $sd) {
            $min = round($median - 2 * $sd, $this->precision);

            foreach ($x as $value) {
                if ($value < $min) {
                    $outliers[] = $value;
                }
            }
        }

        if ($max > $median + 2 * $sd) {
            $max = round($median + 2 * $sd, $this->precision);

            foreach ($x as $value) {
                if ($value > $max) {
                    $outliers[] = $value;
                }
            }
        }

        return array('min' => $min, 'q1' => $q1, 'median' => $median, 'q3' => $q3, 'max' => $max, 'outliers' => $outliers);
    }

    public function hist($x, $bins = false)
    {
        if ($bins === false) {
            $bins = ceil(sqrt(count($x)));
        }

        $max = max($x);
        $min = min($x);
        $bar = ($max - $min) / $bins;

        $bars   = array();
        $ranges = array();

        for ($i = 0; $i < $bins; $i++) {
            $bars[$i] = 0;

            $low  = round($min + $i * $bar, $this->precision);
            $high = round($min + ($i + 1) * $bar, $this->precision);

            $ranges[$i] = "$low-$high";
        }

        foreach ($x as $value) {
            if ($value == $max) {
                $bars[$bins - 1]++;
            } else {
                $bars[floor(($value - $min) / $bar)]++;
            }
        }

        return (array_combine($ranges, $bars));
    }

    public function qqnorm($y)
    {
        asort($y);

        $qqx = array();
        $qqy = array();

        $n = count($y);
        $i = 1;

        foreach ($y as $k => $v) {
            $z = $this->inverseNormCDF(($i - 0.5) / $n);
            $qqx[$k] = $z;
            $qqy[$k] = $v;
            $i++;
        }

        ksort($qqx);
        ksort($qqy);

        $qq = array();

        $qq['x'] = array_values($qqx);
        $qq['y'] = array_values($qqy);

        return $qq;
    }

    /**
     * A ternary plot is a barycentric plot on three variables which sum to a constant. 
     * It graphically depicts the ratios of the three variables as positions in an equilateral 
     * triangle. It is used to show the compositions of systems composed of three species.
     *
     * @param array $x List of float values corresponding to the X sample of a population.
     * @param array $y List of float values corresponding to the Y sample of a population.
     * @param array $z List of float values corresponding to the Z sample of a population.
     *                    
     * @return array Returns x,y coordinates for each data point to render requested ternary plot.
     * @author Khaled Al-Sham'aa <khaled.alshamaa@gmail.com>
     * @url http://en.wikipedia.org/wiki/Ternary_plot
     */
    public function ternary($x, $y, $z)
    {
        $n = count($x);

        $x2d = array();
        $y2d = array();

        for ($i = 0; $i < $n; $i++) {
            $x2d[] = round(0.5 * (2 * $y[$i] + $z[$i]) / ($x[$i] + $y[$i] + $z[$i]), $this->precision);
            $y2d[] = round((sqrt(3) / 2) * $z[$i] / ($x[$i] + $y[$i] + $z[$i]), $this->precision);
        }

        $tri = array();

        $tri['x'] = $x2d;
        $tri['y'] = $y2d;

        return $tri;
    }

    /**
     * Returns the inverse of the standard normal cumulative distribution. 
     * The distribution has a mean of zero and a standard deviation of one.
     * This is an implementation of the algorithm published at:
     * http://home.online.no/~pjacklam/notes/invnorm/
     *
     * @param float $p Is a probability corresponding to the normal distribution between 0 and 1.
     *
     * @return float Inverse of the standard normal cumulative distribution, with a probability of $p
     * @author Khaled Al-Sham'aa <khaled.alshamaa@gmail.com>
     */
    public function inverseNormCDF($p)
    {
        /* coefficients for the rational approximants for the normal probit: */
        $a1    = -3.969683028665376e+01;
        $a2    =  2.209460984245205e+02;
        $a3    = -2.759285104469687e+02;
        $a4    =  1.383577518672690e+02;
        $a5    = -3.066479806614716e+01;
        $a6    =  2.506628277459239e+00;
        $b1    = -5.447609879822406e+01;
        $b2    =  1.615858368580409e+02;
        $b3    = -1.556989798598866e+02;
        $b4    =  6.680131188771972e+01;
        $b5    = -1.328068155288572e+01;
        $c1    = -7.784894002430293e-03;
        $c2    = -3.223964580411365e-01;
        $c3    = -2.400758277161838e+00;
        $c4    = -2.549732539343734e+00;
        $c5    =  4.374664141464968e+00;
        $c6    =  2.938163982698783e+00;
        $d1    =  7.784695709041462e-03;
        $d2    =  3.224671290700398e-01;
        $d3    =  2.445134137142996e+00;
        $d4    =  3.754408661907416e+00;

        $p_low  = 0.02425;
        $p_high    = 1.0 - $p_low;

        if (0 < $p && $p < $p_low) {
            /* rational approximation for the lower region */
            $q = sqrt(-2 * log($p));
            $x = ((((($c1 * $q + $c2) * $q + $c3) * $q + $c4) * $q + $c5) * $q + $c6) / (((($d1 * $q + $d2) * $q + $d3) * $q + $d4) * $q + 1);
        } elseif ($p_low <= $p && $p <= $p_high) {
            /* rational approximation for the central region */
            $q = $p - 0.5;
            $r = $q * $q;
            $x = ((((($a1 * $r + $a2) * $r + $a3) * $r + $a4) * $r + $a5) * $r + $a6) * $q / ((((($b1 * $r + $b2) * $r + $b3) * $r + $b4) * $r + $b5) * $r + 1);
        } else {
            /* rational approximation for the upper region */
            $q = sqrt(-2 * log(1 - $p));
            $x = - ((((($c1 * $q + $c2) * $q + $c3) * $q + $c4) * $q + $c5) * $q + $c6) / (((($d1 * $q + $d2) * $q + $d3) * $q + $d4) * $q + 1);
        }

        return $x;
    }
}
